# -*- coding: utf-8 -*-
"""
Created on Tue Jun  6 11:49:30 2017

@author: Leon
"""

# -*- coding: utf-8 -*-
"""This module is for loading binary image file."""
import cv2
import numpy as np
import struct


def readBsq(imgName, numRow, numCol, numBand, dataFormat):
    """Read an binary image with BSQ format."""
    inputImg = open(imgName, 'rb')
    bytePerPix = struct.calcsize(dataFormat)    # Get the size of data format

    outputImg = np.empty((numRow, numCol, numBand))
    for band in range(numBand):
        for row in range(numRow):
            for col in range(numCol):
                #print(inputImg.read(bytePerPix))
                data = inputImg.read(bytePerPix)    # Read pixel binary value
                
                try:
                #  Unpack and store the binary value to output image array
                    outputImg[row, col, band] = struct.unpack(dataFormat, data)[0]
                except:
                    pass
                
    return outputImg


def main():
    # Required image information
    inputImgName = 'img/MS.img'
    outputImgName = 'img/MS.tif'
    numRow = 800
    numCol = 1200
    numBand = 4
    dataFormat = 'H'        # 2-byte unsigned integer

    outputImg = readBsq(inputImgName, numRow, numCol, numBand, dataFormat)
    for i in range(numBand):
        # Determine the target value range after enhancement
        hist, _ = np.histogram(outputImg[:, :, i], bins=4000)
        accu = hist.cumsum()
        pc = 100. * accu / hist.sum()   # Cumulative probability array
        idxMin = (np.abs(pc - 2)).argmin()      # Index of %2
        idxMax = (np.abs(pc - 98)).argmin()     # Index of %98

        # Max/min value of output image
        valMin = np.where(accu == accu[idxMin])[0][-1]
        valMax = np.where(accu == accu[idxMax])[0][0]

        # Compute parameters of linear transform function
        slope = 255.0 / (valMax - valMin)
        incpt = 0 - valMin * slope

        # Update band values
        outputImg[:, :, i] = outputImg[:, :, i] * slope + incpt

    # Limit the output image values
    outputImg = np.clip(outputImg, 0, 255)
    B = outputImg[:,:,0]/1.5
    G = outputImg[:,:,1]
    R = outputImg[:,:,2]
    in_R = outputImg[:,:,3]
    
    img = np.dstack((B,G,R))
    print(img.shape)
    # Show the result
    
    cv2.imshow(outputImgName, img.astype(np.uint8))
    cv2.waitKey(0)
    cv2.destroyWindow(outputImgName)

    # Save image
    cv2.imwrite(outputImgName, img.astype(np.uint8))


if __name__ == '__main__':
    main()